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ABSTRACT 

In biomolecular systems (especially all-atom models) with many degrees of freedom 
such as proteins and nucleic acids, there exist an astronomically large number of local- 
minimum-energy states. Conventional simulations in the canonical ensemble are of little 
use, because they tend to get trapped in states of these energy local minima. Enhanced 
conformational sampling techniques are thus in great demand. A simulation in general- 
ized ensemble performs a random walk in potential energy space and can overcome this 
difficulty. From only one simulation run, one can obtain canonical-ensemble averages of 
physical quantities as functions of temperature by the single-histogram and/or multiple- 
histogram reweighting techniques. In this article we review uses of the generalized- 
ensemble algorithms in biomolecular systems. Three well-known methods, namely, mul- 
ticanonical algorithm, simulated tempering, and replica-exchange method, are described 
first. Both Monte Carlo and molecular dynamics versions of the algorithms are given. 
We then present various extensions of these three generalized-ensemble algorithms. The 
effectiveness of the methods is tested with short peptide and protein systems. 
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1 INTRODUCTION 



Conventional Monte Carlo (MC) and molecular dynamics (MD) simulations of biomolecules 
are greatly hampered by the multiple-minima problem. The canonical fixed-temperature 
simulations at low temperatures tend to get trapped in a few of a huge number of local- 
minimum-energy states which are separated by high energy barriers. One way to overcome 
this multiple-minima problem is to perform a simulated annealing (SA) simulation pQ, 
and it has been widely used in biomolecular systems (see, e.g., Refs. [2]-[5] for earlier 
applications). The SA simulation mimicks the crystal- making process, and temperature 
is lowered very slowly from a sufficiently high temperature to a low one during the SA 
simulation. The Boltzmann weight factor is dynamically changed, and so the thermal 
equilibrium is continuously broken. Hence, although the global-minimum potential en- 
ergy or the value close to it may be found, accurate thermodynamic averages for fixed 
temperatures cannot be obtained. 

A class of simulation methods, which are referred to as the generalized- ensemble algo- 
rithms, overcome both above difficulties, namely the multipole-minima problem and inac- 
curate thermodynamic averages (for reviews see, e.g., Refs. In the generalized- 
ensemble algorithm, each state is weighted by an artificial, non-Boltzmann probability 
weight factor so that a random walk in potential energy space may be realized. The 
random walk allows the simulation to escape from any energy barrier and to sample much 
wider conformational space than by conventional methods. Unlike SA simulations, the 
weight factors are fixed during the simulations so that the eventual reach to the thermal 
equilibrium is guaranteed. From a single simulation run, one can obtain accurate en- 
semble averages as functions of temperature by the single-histogram [T7] and/or multiple- 
histogram [m [19] reweighting techniques (an extension of the multiple-histogram method 
is also referred to as the weighted histogram analysis method (WHAM) p!9]). 

One of the most well-known generalized-ensemble algorithms is perhaps the multi- 
canonical algorithm (MUCA) [201121] (for reviews see, e.g., Refs. |221I23])- The method 
is also referred to as entropic sampling [25] and adaptive umbrella sampling [SB] of 
the potential energy [27]. MUCA can also be considered as a sophisticated, ideal realiza- 
tion of a class of algorithms called umbrella sampling [28]. Also closely related methods 
are Wang-Landau method [291 [20], which is also referred to as density of states Monte 
Carlo [21], and Meta Dynamics [32] • See also Ref. [33] • MUCA and its generalizations 
have been applied to spin systems (see, e.g., Refs. [31]^[in])- MUCA was also introduced 
to the molecular simulation field [H]. Since then MUCA and its generalizations have 
been extensively used in many applications in protein and other biomolecular systems 
[l2]-[76]. Molecular dynamics version of MUCA has also been developed [HI [521 [27] 
(see also Refs. [771 [IH] for the Langevin dynamics version). MUCA has been extended 
so that fiat distributions in other variables instead of potential energy may be obtained 
(see, e.g., Refs. [351 [Ml [13 [531 [SSI [HHl [H]). This can be considered as a special case of 
the multidimensional (or, multivariable) extensions of MUCA, where a multidimensional 
random walk in potential energy space and in other variable space is realized (see, e.g., 
Refs. [m [S31 [SH [701 [TB] ) • In this article, we just present one of such methods, namely, the 
multibaric-multithermal algorithm (MUBATH) where a two-dimensional random walk in 
both potential energy space and volume space is realized [7D]-[75]. 

While a simulation in multicanonical ensemble performs a free ID random walk in 
potential energy space, that in simulated tempering (ST) [751 [ZH] (the method is also 
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referred to as the method of expanded ensemble [78]) performs a free random walk in 
temperature space (for a review, see, e.g., Ref. [50] )• This random walk, in turn, induces 
a random walk in potential energy space and allows the simulation to escape from states 
of energy local minima. ST and its generalizations have also been applied to chemical 
physics field and biomolecular systems [8ll[82lll9l[50l|83l[8l[85l|86l[HIll8i 

MUCA and ST are powerful, but the probability weight factors are not a priori known 
and have to be determined by iterations of short trial simulations. This process can be 
non-trivial and very tedius for complex systems with many degreees of freedom. 

In the replica- exchange method (REM) [Hn]^[H2], the difficulty of weight factor deter- 
mination is greatly alleviated. (A closely related method was independently developed in 
Ref. [93]. Similar methods in which the same equations are used but emphasis is laid on 
optimizations have been developed [M| [95]. REM is also referred to as multiple Markov 
chain method [HS] and parallel tempering [5U]. Details of literature about REM and re- 
lated algorithms can be found in recent reviews [TU| [T3].) In this method, a number 
of non-interacting copies (or, replicas) of the original system at different temperatures 
are simulated independently and simultaneously by the conventional MC or MD method. 
Every few steps, pairs of replicas are exchanged with a specified transition probability. 
The weight factor is just the product of Boltzmann factors, and so it is essentially known. 

REM has already been used in many applications in protein systems |98]- |113] . Other 
molecular simulation fields have also been studied by this method in various ensembles 
|114] - [TT8] . Moreover, REM was introduced to the quantum chemistry field |119] . The 
details of molecular dynamics algorithm for REM, which is referred to as the Replica- 
Exchange Molecular Dynamics (REMD) have been worked out in Ref. [HI], and this led 
to a wide application of REM in the protein folding and related problems (see, e.g., 
Refs. [I20]-[I12]). 

However, REM also has a computational difficulty: As the number of degrees of 
freedom of the system increases, the required number of replicas also greatly increases, 
whereas only a single replica is simulated in MUCA and ST. This demands a lot of com- 
puter power for complex systems. Our solution to this problem is: Use REM for the weight 
factor determinations of MUCA, which is much simpler than previous iterative methods 
of weight determinations, and then perform a long MUCA production run. The method is 
referred to as the replica- exchange multicanonical algorithm (REMUCA) |103[ I107[ 1108] . 
In REMUCA, a short replica-exchange simulation is performed, and the multicanonical 
weight factor is determined by the multiple-histogram reweighting techniques [TSl [TU] . 
Another example of a combination of REM and ST is the replica- exchange simulated tem- 
pering (REST) [SI]. In REST, a short replica-exchange simulation is performed, and the 
simulated tempering weight factor is determined by the multiple-histogram reweighting 
techniques [TH [19] . 

We have introduced two further extensions of REM, which we refer to as multicanonical 
replica-exchange method (MUCAREM) [1031 [Ml [IDS] (see also Refs. [Ml [M]) and 
simulated tempering replica- exchange method (STREM) [85] (see also Ref. |133] for a 
similar idea). In MUCAREM, a replica-exchange simulation is performed with a small 
number of replicas each in multicanonical ensemble of different energy ranges. In STREM, 
on the other hand, a replica-exchange simulation is performed with a small number of 
replicas in "simulated tempering" ensemble of different temperature ranges. 

Finally, one is naturally led to a multidimensional (or, multivariable) extension of 
REM, which we refer to as the multidimensional replica- exhcange method (MREM) |101] 
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(see also Refs. |143[ 1115] ). (The method is also referred to as generalized parallel sampling 
[144] . Hamiltonian replica- exchange method |1U6] . and Model Hoppina \\'ib\ .) Some other 
special cases of MREM can be found in, e.g., Refs. [M ESDI EZl [M [M [M fTH] . 
Another special realization of MREM is replica- exchange umbrella sampling (REUS) |101j 
and it is particularly useful in free energy calculations (see also Ref. |102] for a similar 
idea). In this article, we just present one of such methods, namely, the replica-exchange 
method in the isobaric-isothermal ensemble, where not only temperature values but also 
pressure values are exchanged in the replica-exchange processes |116l I118[ [TT| 11261 1127] . 
(The results of the first such application of the two-dimensional replica-exchange simula- 
tions in the isobaric-isothermal ensemble were presented in Ref. [TT].) 

In this article, we describe the generalized-ensemble algorithms mentioned above. 
Namely, we first review the three familiar methods: REM, ST, and MUCA. We then 
describe various extensions of these methods |1U11 1152|, 11531 11541 1155j . Examples of the 
results by some of these algorithms are then presented. 



2 GENERALIZED-ENSEMBLE ALGORITHMS 

2.1 Replica- Exchange Method 

Let us consider a system of N atoms of mass {k = 1, ■ ■ ■ , A^) with their coordinate 
vectors and momentum vectors denoted by g = {qi, ■ ■ ■ ,q^} and p = {p^, ■ ■ ■ ,pj^}, 
respectively. The Hamiltonian H{q,p) of the system is the sum of the kinetic energy 
K{p) and the potential energy E{q): 

H{q,p)=K{p) + E{q) , (1) 

where 

N 2 

K{p) = E ^ • (2) 

In the canonical ensemble at temperature T each state x = {q,p) with the Hamiltonian 
H[q,p) is weighted by the Boltzmann factor: 

WB{x-T)=exp{-f3H{q,p)) , (3) 

where the inverse temperature /3 is defined by /3 = 1/kBT {kB is the Boltzmann constant). 
The average kinetic energy at temperature T is then given by 

Because the coordinates q and momenta p are decoupled in Eq. ([T]), we can suppress 
the kinetic energy part and can write the Boltzmann factor as 

Wnix; T) = Wb{E; T) = exp(-/3E) . (5) 

The canonical probability distribution of potential energy Pnvt(-E'; T) is then given by 
the product of the density of states n{E) and the Boltzmann weight factor Wb{E\ T): 

Pnvt(^; T) oc n{E)WB[E- T) . (6) 
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Because n{E) is a rapidly increasing function and the Boltzmann factor decreases ex- 
ponentially, the canonical ensemble yields a bell-shaped distribution of potential energy 
which has a maximum around the average energy at temperature T. The conventional 
MC or MD simulations at constant temperature are expected to yield Pnvt(-E'; 7"). A 
MC simulation based on the Metropolis algorithm |156] is performed with the following 
transition probability from a state x of potential energy E to a. state x' of potential energy 
E': 

w{x ^ x') = min (^1, = min (1, exp (-/3Ai?)) . (7) 

where 

AE = E' -E . (8) 

A MD simulation, on the other hand, is based on the following Newton equations of 
motion: 

9. = — , (9) 

mfc 

OE 

P. = -g^ = A.. (10) 

where ff, is the force acting on the k-th atom {k = 1,---,N). This set of equations 
actually yield the microcanonical ensemble, however, and we have to add a thermostat 
in order to obtain the canonical ensemble at temperature T. Here, we just follow Nose's 
prescription |157[ 1158] , and we have 

= ^, (11) 
rrik 

d^Ej s s 
i = «f . (13) 

A'' 2 

Ps = E — - ^N^bT = 3NkB (T(t) - T) , (14) 



k-- 



where s is Nose's scaling parameter, P, is its conjugate momentum, Q is its mass, and 
the "instantaneous temperature" T(t) is defined by 

(15) 

However, in practice, it is very difficult to obtain accurate canonical distributions of 
complex systems at low temperatures by conventional MC or MD simulation methods. 
This is because simulations at low temperatures tend to get trapped in one or a few of 
local-minimum-energy states. This difficulty is overcome by, for instance, the generalized- 
ensemble algorithms, which greatly enhance conformational sampling. 

The replica- exchange method (REM) is one of effective generalized-ensemble algo- 
rithms. The system for REM consists of M non-interacting copies (or, replicas) of the 
original system in the canonical ensemble at M different temperatures (m = 1, ■ ■ ■ , M). 
We arrange the replicas so that there is always exactly one replica at each temperature. 
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Then there exists a one-to-one correspondence between rephcas and temperatures; the 
label i [i = 1, ■ ■ ■ , M) for replicas is a permutation of the label m [m = 1, ■ ■ ■ , M) for 
temperatures, and vice versa: 

1 m = m(i) = f~^{i) , 

where f{m) is a permutation function of m and f~^{i) is its inverse. 

Let X = ^Xi^^^\ ■ ■ ■ , = ■ ■ ■ , a:[^jy-^| stand for a "state" in this general- 

ized ensemble. Each "substate" x]^ is specified by the coordinates g'*' and momenta p'*' 
of atoms in replica i at temperature T^: 

xW^(gW,p[^l)^ . (17) 

Because the replicas are non-interacting, the weight factor for the state X in this 
generalized ensemble is given by the product of Boltzmann factors for each replica (or at 
each temperature): 



M M 



i=l m=l 



:i8) 



where i{m) and m(i) are the permutation functions in Eq. (ITB]) . 

We now consider exchanging a pair of replicas in this ensemble. Suppose we exchange 
replicas i and j which are at temperatures and T„, respectively: 

X — ^ - ■ ■ , , ■ ■ ■ , ' ' '} ^ ^ ~ {' ' ' ' '^'^^ ' ' ' ' ' -^ri ''''}■ i^^) 

Here, i, j, m, and n are related by the permutation functions in Eq. fll6l) . and the exchange 
of replicas introduces a new permutation function /': 

/ i = fi^) >J= f'i^) , (nn\ 

1 J = f{n) I = f\n) . ^^^^ 



The exchange of replicas can be written in more detail as 



V / n V / n 



(21) 



where the definitions for pt*^' and p^^^' will be given below. We remark that this process is 
equivalent to exchanging a pair of temperatures and T„ for the corresponding replicas 
i and j as follows: 



\ / n \ /I 



(22) 



In the original implementation of the replica- exchange method (REM) [90]-[92], Monte 
Carlo algorithm was used, and only the coordinates q (and the potential energy function 
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E[q)) had to be taken into account. In molecular dynamics algorithm, on the other 
hand, we also have to deal with the momenta p. We proposed the following momentum 
assignment in Eq. (ETj) (and in Eq. 0221) ) 



pi' 




V 



{]]' 



p. 



Tm \j] 

T 

■J- n. 




(23) 



which we believe is the simplest and the most natural. This assignment means that we 
just rescale uniformly the velocities of all the atoms in the replicas by the square root of 
the ratio of the two temperatures so that the temperature condition in Eq. may be 
satisfied immediately after replica exchange is accepted. 

The transition probability of this replica-exchange process is given by the usual Metropo- 
lis criterion: 



w{X X') 



W X, 



X, 



mm 



W"rem(^0 



min (1, exp (—A)) 



(24) 



where in the second expression (i.e., w{x^^\x^^)) we explicitly wrote the pair of replicas 
(and temperatures) to be exchanged. From Eqs. ([I]), ([2]), flTH]) . and f l2^ . we have 



W^rem(XO 
H^rem(X) 



exp [-P^ [k (pbl') + E (gbl)] - /3„ [k (pH') + E {q 

+(3^^ (pW) + E (q^A] + (3n [k (pb-l) + E (gb-])] } , 



(25) 



(3^ \e (gbl) -E(q 



Because the kinetic energy terms in this equation all cancel out, A in Eq. fl24|l becomes 



= {Pm-Pn)(E(q^^^)-E(q 



(26) 
(27) 



Here, j, m, and n are related by the permutation functions in Eq. ( 1T6|) before the replica 
exchange: 

\j =f{n). ^28) 

Note that after introducing the momentum rescaling in Eq. f l23|) . we have the same 
Metropolis criterion for replica exchanges, i.e., Eqs. flMj) and fl27|) . for both MC and MD 
versions. 

Without loss of generality we can assume Ti < T2 < ■ ■ ■ < Tm- The lowest temperature 
Ti should be sufficiently low so that the simulation can explore the global-minimum-energy 
region, and the highest temperature Tm should be sufficiently high so that no trapping 
in an energy-local- minimum state occurs. A simulation of the replica- exchange method 
(REM) is then realized by alternately performing the following two steps: 

1. Each replica in canonical ensemble of the fixed temperature is simulated simultaneously 
and independently for a certain MC or MD steps. 
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2. A pair of replicas at neighboring temperatures, say and are exchanged 



with the probabihty w [x^^ ^^m+i ) ™- Eq- 

Note that in Step 2 we exchange only pairs of replicas corresponding to neighboring tem- 
peratures, because the acceptance ratio of the exchange process decreases exponentially 
with the difference of the two /3's (see Eqs. fl27p and f l24p ) . Note also that whenever a 
replica exchange is accepted in Step 2, the permutation functions in Eq. f ll6p are updated. 
A random walk in "temperature space" is realized for each replica, which in turn induces 
a random walk in potential energy space. This alleviates the problem of getting trapped 
in states of energy local minima. 

The REM simulation is particularly suitable for parallel computers. Because one can 
minimize the amount of information exchanged among nodes, it is best to assign each 
replica to each node (exchanging pairs of temperature values among nodes is much faster 
than exchanging coordinates and momenta). This means that we keep track of the per- 
mutation function m{i\ t) = f~^{i] t) in Eq. f|T6|) as a function of MC or MD step t during 
the simulation. After parallel canonical MC or MD simulations for a certain steps (Step 
1), M/2 pairs of replicas corresponding to neighboring temperatures are simulateneously 
exchanged (Step 2), and the pairing is alternated between the two possible choices, i.e., 
(Ti,T2), (T3,T4), ■■■ and {T^.T^), (T.^T,),---. 

After a long production run of a replica-exchange simulation, the canonical expectation 
value of a physical quantity A at temperature {m = 1, ■ ■ ■ , M) can be calculated by 
the usual arithmetic mean: 

1 

<A>T^=—Y.A{xUk)) , (29) 

where Xm{k) {k = 1, - ■ ■ the configurations obtained at temperature and 

is the total number of measurements made at T = T^. The expectation value at any 
intermediate temperature T (= l/ksf^) can also be obtained as follows: 

J2 A(E)Pnvt(^;T) A{E)n{E)exp{-f3E) 

< A >T= = . (30) 

P^yT{E;T) Y n{E)exp{-l3E) 

E E 

The summation instead of integration is used in Eq. fl5U]l . because we often discretize the 
potential energy E with step size e [E = Ei\i = 1,2,- ■ ■). Here, the explicit form of the 
physical quantity A should be known as a function of potential energy E. For instance, 
A{E) = E gives the average potential energy < E >t as a function of temperature, and 
A{E) = 13'^ {E- < E >tY gives specific heat. 

The density of states n{E) in Eq. (!30|) is given by the multiple-histogram reweighting 
techniques [ISl US] as follows. Let N^i^E) and be respectively the potential-energy 
histogram and the total number of samples obtained at temperature = 1/ ^B/^m (^ = 
1, ■ ■ ■ , M). The best estimate of the density of states is then given by [T8| [19] 

M 

n{E) = , (31) 

Y 9^ exp(/„ - I3^E) 

m=l 
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where we have for each m (= 1, ■ ■ ■ , M) 

exp(-/„) = ^ n{E) exp{-(3mE) . (32) 

E 

Here, Qm = ^ + Sr^, and is the integrated autocorrelation time at temperature T^. 
For many systems the quantity Qm can safely be set to be a constant in the reweighting 
formulae [19], and hereafter we set Qm = ^■ 

Note that Eqs. ( 13T|) and ( 132|) are solved self-consistently by iteration [HI [19] to obtain 
the density of states n{E) and the dimensionless Helmholtz free energy f^- Namely, we 
can set all the (m = 1, ■ ■ ■ , M) to, e.g., zero initially. We then use Eq. (EI]) to obtain 
n{E), which is substituted into Eq. (132|1 to obtain next values of /„, and so on. 

Moreover, the ensemble averages of any physical quantity A (including those that 
cannot be expressed as functions of potential energy) at any temperature T (= l/ksP) can 
now be obtained from the "trajectory" of configurations of the production run. Namely, 
we first obtain fm [m = 1, ■ ■ ■ , M) by solving Eqs. flHT]) and fl5^ self-consistently, and 
then we have |1U7] (see also |159] ) 

M rim 1 

Y: Y.M^m{k))- exp[-PE{xUm 

< ^ >-= ir^m 1 ' (33) 

E E -I ^^p [-(3E{xm{m 

m=i k=i j2 rif, exp [f(, - (3(,E{xmik))] 

where Xm{k) {k = 1, ■ ■ ■ , rim) are the configurations obtained at temperature T^,. 

Eqs. (13 op and (13 ip or any other equations which involve summations of exponential 
functions often encounter with numerical difficulties such as overflows. These can be 
overcome by using, for instance, the following equation [221 HSU]: For C = A + B (with 
A > and B > 0) we have 



InC =ln 



, , , / mini A, B) 
max(A,5) 1 ' ^ ' ' 



max(A, B) 



(34) 



max(ln/l. In 5) + In {1 + exp [min(ln A, In B) — max(ln A, In 5)]} . 



We now give more details about the momentum rescaling in Eq. fl23|) |161] . Actually, 
Eq. (123]) is only valid for the Langevin dynamics |162] . Andersen thermostat [163] . and 
Gaussian constraint method [1641 11651 1166] . The former two thermostats are based on 
the weight factor of Eq. ^ with Eqs. ([T]) and ([2]), and the Gaussian contraint method is 
based on the following weight factor: 

W{q.p) = 6(j:^- ^) exp [-mq)] ■ (35) 



\k-. 



For Nose's method [1571 1158] . which gives the equations of motion in Eqs. (fTT|) - ([ni) . 
the Nose Hamiltonian is given by 

^Nose = E + E{q) + ^ + gk^Tlog s. (36) 
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Here, g (= 3iV) is the number of degrees of freedom, s is a position variable of the 
thermostat, Pg is a momentum conjugate to s, and is a virtual momentum, which is 
related to the real momenta p^, as = p^,,/ s. The weight factor for the Nose's method is 
then given by 

W{q,p,s,Ps) = 5{H^osc-£) , (37) 

where £ is the initial value of H-^ose- Namely, in the Nose's method, the entire system 
including the thermostat is in the microcanonical ensemble. Note that the mass Q of the 
thermostat can have different values in each replica in REMD simulations. The rescaling 
method for the Nose thermostat is given by Eq. fl2^ and 



' ^nQn 
TmQm 



p{i] p[0]' 
s 1 s 



J exp 



s'-'l exp 



1 (E{q^^)-£r, 



I TmQm 

E{ql^ 



.plJ] 



T 



1 /E(gl^J)-^„ E{q^^^)-£, 



T 



T 

J. rr 



(3J 



(39) 



where £m and £„, are the initial values of -f^Nose in the simulations with and T„, 
respectively, before the replica exchange. Note that the real momenta have to be used in 
the rescaling method in Eq. (125]) . not the virtual momenta. 

For the Nose-Hoover thermostat |167j . the states are specified by the following weight 
factor: 



^"(^,^,0 = exp 



(40) 



where C is a velocity of the thermostat and Q is its mass parameter. The rescaling method 
for the Nose-Hoover thermostat is given by Eq. fl2^ and 



TmQn 



(41) 



where Qm and Q„ are the mass parameters in the replicas at temperature values and 
T„, respectively, before the replica exchange. 

The rescaling method for the Nose-Hoover thermostat can be generalized to the Nose- 
Hoover chains |168] in a similar way. The weight factor for the Nose-Hoover chains is 
given by 



W{q,pXi 



exp 



£=1 ^ /. 



(42) 



1, ■ ■ ■ , £) is the velocity of the £-th ther- 



where C is the number of thermostats, Q {i 
mostat, and Qi (£= 1, ■■■,£) is a mass parameter corresponding to the i-th thermostat. 
A rescaling method for REMD with the Nose-Hoover chains is given by Eq. (l23l) and the 
following: 



r 



Cf, (^ 



X), 



(43) 



where Qm,i and Qn,e. are the mass parameters in the replicas at temperature values T, 
and T„, respectively, which correspond to the £-th thermostat. 
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In the Nose-Poincare thermostat [169] . the states are specified by x = (g, i?, s, Pg) and 
the weight factor is given by 

W{q, p, s, Ps) oc 6 [s (i/Noso - S)] , (44) 

where i^Nosc is the Nose Hamiltonian in Eq. fl36|) and S is its initial value. A rescaling 
method of the Nose-Poincare thermostat is the same as in the Nose's thermostat and 
given by Eqs. 022]), dM]), and (EH]) above. 



2.2 Simulated Tempering 

We now introduce another generalized-ensemble algorithm, the simulated tempering (ST) 
method [781 EHj. In this method temperature itself becomes a dynamical variable, and 
both the configuration and the temperature are updated during the simulation with a 
weight: 

Wst{E; T) = exp (-/3E + a(T)) , (45) 

where the function a{T) is chosen so that the probability distribution of temperature is 
fiat: 

PgT(T) = j dE n{E) WsTiE; T) = J dE n{E) exp {-(3E + a(T)) = constant . (46) 

Hence, in simulated tempering, temperature is sampled uniformly. A free random walk in 
temperature space is realized, which in turn induces a random walk in potential energy 
space and allows the simulation to escape from states of energy local minima. 

In the numerical work we discretize the temperature in M different values, (m = 
1, • • • , M). Without loss of generality we can order the temperature so that Ti < T2 < 
■ ■ ■ < Tm- The probability weight factor in Eq. ( H5|l is now written as 



Wst{E- = exp(-/3^E + a J , (47) 
where = ci{Tm) (m = 1, ■ ■ ■ , M). Note that from Eqs. fllB]) and f HT]) we have 

exp(-am) (X J dE n{E) exp(-/3mE) . (48) 

The parameters are therefore "dimensionless" Helmholtz free energy at temperature 
Tm (i.e., the inverse temperature (3m multiplied by the Helmholtz free energy). 

Once the parameters determined and the initial configuration and the initial 

temperature Tm are chosen, a simulated tempering simulation is realized by alternately 
performing the following two steps [7S1 US\ '■ 

1. A canonical MC or MD simulation at the fixed temperature Tm is carried out for a 
certain steps. 

2. The temperature Tm is updated to the neighboring values Tm±i with the configura- 
tion fixed. The transition probability of this temperature- up dating process is given 
by the following Metropolis criterion (see Eq. (H7|) ): 

w{Tm ^ T„,±0 = min (^1, ^^^0) = (1^ ^xp (- A)) , (49) 

where 

A = i(3m±i - (3m) E - {am±i - am) ■ (50) 
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Note that in Step 2 we update the temperature only to the neighboring temperatures in 
order to secure sufficiently large acceptance ratio of temperature updates. 

We remark that when MD simulations are performed in Step 1 above, we also have to 
deal with the momenta p, and the kinetic energy term should be included in the weight 
factor. When temperature Tmo±i is accepted for T-update in Step 2, we rescale the 
momenta in the same way as in REMD [99l 11531 1155j : 



The kinetic energy terms then cancel out in Eq. (l50l) and we can use the same A in 
the Metropolis criterion in Step 2 for both MC and MD simulations. Similar momentum 
scaling formulae given above should also be introduced for various other thermostats [161] . 

The simulated tempering parameters = ci{Tm) (m = 1, ■ ■ ■ , M) can be determined 
by iterations of short trial simulations (see, e.g., Refs. [80l[82l|50] for details). This process 
can be non-trivial and very tedius for complex systems. 

After the optimal simulated tempering weight factor is determined, one performs a long 
simulated tempering run once. The canonical expectation value of a physical quantity A 
at temperature [m = 1, ■ ■ ■ , M) can be calculated by the usual arithmetic mean from 
Eq. (l29l) . The expectation value at any intermediate temperature can also be obtained 
from Eq. (l30l) . where the density of states is given by the multiple-histogram reweighting 
techniques [181 US]- Namely, let Nm{E) and Um be respectively the potential-energy 
histogram and the total number of samples obtained at temperature = 1/ kB/3m ijn = 
1, ■ ■ ■ , M). The best estimate of the density of states is then given by solving Eqs. (PT]) 
and fl32|l self-consistently. 

Moreover, the ensemble averages of any physical quantity A (including those that 
cannot be expressed as functions of potential energy) at any temperature T (= I/Zcb/?) 
can now be obtained from Eq. fl33|) . 

2.3 Multicanonical Algorithm 

The third generalized-ensemble algorithm that we present is the multicanonical algo- 
rithm (MUCA) [201 EI]- In the multicanonical ensemble, each state is weighted by a 
non-Boltzmann weight factor Wmvca{E) (which we refer to as the multicanonical weight 
factor) so that a uniform potential energy distribution Pmuca{E) is obtained: 



The flat distribution implies that a free one-dimensional random walk in the potential 
energy space is realized in this ensemble. This allows the simulation to escape from any 
local minimum-energy states and to sample the configurational space much more widely 
than the conventional canonical MC or MD methods. 

The definition in Eq. flS^ implies that the multicanonical weight factor is inversely 
proportional to the density of states, and we can write it as follows: 




(51) 



^muca(^) oc n{E)WMucA{E) = constant . 



(52) 



Wmvca{E) = exp [-(3aEMucA{E; T,)] 



1 



(53) 



n{E) 
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where we have chosen an arbitrary reference temperature, = l/k-Qf^a, and the ^^multi- 
canonical potential energt/' is defined by 



EMVCA{E;Ta) = kBTJnn{E) = TaS{E) . (54) 

Here, S{E) is the entropy in the microcanonical ensemble. Because the density of states 
of the system is usually unknown, the multicanonical weight factor has to be determined 
numerically by iterations of short preliminary runs f2U\ UT\ . 

A multicanonical MC simulation is performed, for instance, with the usual Metropolis 
criterion |156] : The transition probability of state x with potential energy E to state x' 
with potential energy E' is given by 

w[x x) = mm 1, — — - = mm 1, = mm (1, exp (-/^qAEmuca)) , 



WMUcA{E)j \ 'n{E') 



(55) 



where 



A-Bmuca = EyivcAiE'; Ta) — Emvca{E; Ta) ■ (56) 

The MD algorithm in the multicanonical ensemble also naturally follows from Eq. (1531) . 
in which the regular constant temperature MD simulation (with T = Ta) is performed by 
replacing E by ^muca in Eq. 1^ gHl [52]: 

dEMucA{E;Ta) s (9-Emuca(-E; T^) s 
= d^, = dE (57) 

If the exact multicanonical weight factor WuvcAiE) is known, one can calculate the 
ensemble averages of any physical quantity A at any temperature T (= l/k^P) from 
Eq. (15U]) . where the density of states is given by (see Eq. 



^ ^ Wmvca{E) ^ ^ 

In general, the multicanonical weight factor Wmvca{E) , or the density of states n{E), 
is not a priori known, and one needs its estimator for a numerical simulation. This 
estimator is usually obtained from iterations of short trial multicanonical simulations. 
The details of this process are described, for instance, in Refs. [Ml |13]. However, the 
iterative process can be non-trivial and very tedius for complex systems. 

In practice, it is impossible to obtain the ideal multicanonical weight factor with com- 
pletely uniform potential energy distribution. The question is when to stop the iteration 
for the weight factor determination. Our criterion for a satisfactory weight factor is that 
as long as we do get a random walk in potential energy space, the probability distribu- 
tion Pmuca(-E') does not have to be completely fiat with a tolerance of, say, an order of 
magnitude deviation. In such a case, we usually perform with this weight factor a mul- 
ticanonical simulation with high statistics (production run) in order to get even better 
estimate of the density of states. Let A^muca(-E') be the histogram of potential energy 
distribution Pmuca(-E') obtained by this production run. The best estimate of the density 
of states can then be given by the single-histogram reweighting techniques [17] as follows 
(see the proportionality relation in Eq. (1521) ): 
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By substituting this quantity into Eq. fl30|) . one can calculate ensemble averages of phys- 
ical quantity A{E) as a function of temperature. Moreover, the ensemble averages of any 
physical quantity A (including those that cannot be expressed as functions of potential 
energy) at any temperature T (= l/fce/?) can also be obtained as long as one stores the 
"trajectory" of configurations from the production run. Namely, we have 

ria 

J2 A{xk)W^ljcAiE{xk)) exp [-l3E{xk)] 

< A >T= '-^a ' (60) 

E ^uhcAiEixk)) exp i-m^k)] 

k=l 

where Xk is the configuration at the k-th MC (or MD) step and Ug is the total number of 
configurations stored. Note that when A is a function of E, Eq. ( 160|) reduces to Eq. ( 130|) 
where the density of states is given by Eq. (ISUjl . 

Some remarks are in order. The major advantage of REM over other generalized- 
ensemble methods such as simulated tempering [781 ES] and multicanonical algorithm 
|20[ [2Tj lies in the fact that the weight factor is a priori known (see Eq. f lTSj) ). while in 
simulated tempering and multicanonical algorithm the determination of the weight factors 
can be very tedius and time-consuming. In REM, however, the number of required replicas 
increases greatly (oc \/N) as the system size N increases [50], while only one replica is used 
in simulated tempering and multicanonical algorithm. This demands a lot of computer 
power for complex systems. Moreover, so long as optimal weight factors can be obtained, 
simulated tempering and multicanonical algorithm are more efficient in sampling than the 
replica-exchange method [151 ESI 11081 [75]. 

2.4 Replica-Exchange Simulated Tempering and Replica- Exchange 
Multicanonical Algorithm 

The replica- exchange simulated tempering (REST) [HI] and the replica- exchange multi- 
canonical algorithm (REMUCA) [1031 1107[ 1108] overcome both the difficulties of ST and 
MUCA (the weight factor determinations are non-trivial) and REM (many replicas, or a 
lot of computation time, are required). 

In REST [HI], we first perform a short REM simulation (with M replicas) to deter- 
mine the simulated tempering weight factor and then perform with this weight factor a 
regular ST simulation with high statistics. The first step is accomplished by the multiple- 
histogram reweighting techniques [TS] [19], which give the dimensionless Helmholtz free 
energy. Let Nm{E) and be respectively the potential-energy histogram and the to- 
tal number of samples obtained at temperature (= l/k^f^m) of the REM run. The 
dimensionless Helmholtz free energy fm is then given by solving Eqs. ( 13T]) and ( 132|) self- 
consistently by iteration. 

Once the estimate of the dimensionless Helmholtz free energy fm are obtained, the 
simulated tempering weight factor can be directly determined by using Eq. fl47p where we 
set Qm = fm (compare Eq. (HSj) with Eq. ( 15^ ). A long simulated tempering run is then 
performed with this weight factor. Let Nm{E) and Um be respectively the potential-energy 
histogram and the total number of samples obtained at temperature Tm (= l/fes/Sm) from 
this simulated tempering run. The multiple-histogram reweighting techniques of Eqs. f l3T|) 
and ( |32|) can be used again to obtain the best estimate of the density of states n{E). The 
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expectation value of a physical quantity A at any temperature T (= I/Zcb/?) is then 
calculated from Eq. f l5U]) . 

We now present the replica- exchange multicanonical algorithm (REMUCA) |1U31 llUTl 
llOSj . In REMUCA, just as in REST, we first perform a short REM simulation (with 
M replicas) to determine the multicanonical weight factor and then perform with this 
weight factor a regular multicanonical simulation with high statistics. The first step is 
accomplished by the multiple-histogram reweighting techniques [HI [19] , which give the 
density of states. Let N^i^E) and be respectively the potential-energy histogram and 
the total number of samples obtained at temperature (= l/fce/^m) of the REM run. 
The density of states n{E) is then given by solving Eqs. fl3T|) and fl32l) self-consistently by 
iteration. 

Once the estimate of the density of states is obtained, the multicanonical weight factor 
can be directly determined from Eq. (see also Eq. 0541] ) . Actually, the density of states 
n{E) and the multicanonical potential energy, _E'muca(-E'; To), thus determined are only 
reliable in the following range: 

Ei<E<Em , (61) 



(62) 



where 

< E >T, 

< E >Tm 

and Ti and Tm are respectively the lowest and the highest temperatures used in the REM 
run. Outside this range we extrapolate the multicanonical potential energy linearly |103] : 




<9-E'muca(-E; To) 



^muca(-E') = < 



dE 

Em\jca{E; Tq) , 
dEMvcAiE] To) 



dE 



{E - El) + Emuca(^i; To) , for E < ^i, 

E=Ei 

for Ei<E< Em, 
(E — Em) + -E'muca(-E'm; To) , ioi E > Em- 



(63) 

For Monte Carlo method, the weight factor for REMUCA is given by substituting 
Eq. ([63]) into Eq. ([53]) [11131 [TnTj : 



Wmvca{E) = exp [-f3oSMvcA{E)] 



exp 
1 



-f3iE) , for E < El, 

for Ei<E < Em, 



n{E) ' 

exp {-I3mE) , for E > Em 



(64) 



The multicanonical MC and MD runs are then performed respectively with the Metropo- 
lis criterion of Eq. ([55]) and with the modified Newton equation in Eq. (15 7|) . in which 
^muca(-£') in Eq. ( 163]) is substituted into £'muca(-E'; To). We expect to obtain a flat po- 
tential energy distribution in the range of Eq. (I6T]) . Finally, the results are analyzed by 
the single-histogram reweighting techniques as described in Eq. (159 p (and Eq. ( I30p ). 

The formulations of REST and REMUCA are simple and straightforward, but the 
numerical improvement is great, because the weight factor determination for ST and 
MUCA becomes very difficult by the usual iterative processes for complex systems. 
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2.5 Simulated Tempering Replica-Exchange Method and Mul- 
ticanonical Replica-Exchange Method 

In the previous subsection we presented REST and REMUCA, which use a short REM run 
for the determination of the simulated temepering weight factor and the multicanonical 
weight factor, respectively. Here, we present two modifications of REM and refer to 
the new methods as the simulated tempering replica- exchange method (STREM) [85] and 
multicanonical replica- exchange method (MUCAREM) [M HDZl HOS]. In STREM the 
production run is a REM simulation with a few replicas that perform ST simulations 
with different temperature ranges. Likewise, in MUCAREM the production run is a 
REM simulation with a few replicas in multicanonical ensembles, i.e., different replicas 
perform MUCA simulations with different energy ranges. 

While ST and MUCA simulations are usually based on local updates, a replica- 
exchange process can be considered to be a global update, and global updates enhance 
the conformational sampling further. 

3 MULTIDIMENSIONAL EXTENSIONS OF THE 

THREE GENERALIZED-ENSEMBLE ALGORITHMS 

3.1 General Formulations 

We now give the general formulations for the multidimensional generalized-ensemble al- 
gorithms [1521 1153[ 1154] . Let us consider a generalized potential energy function E^{x), 
which depends on L parameters A = (A^^'', ■ ■ ■ , A^^^), of a system in state x. Although 
E^[x) can be any function of A, we consider the following specific generalized potential 
energy function for simplicity: 

E;^(x) = Eo(x) + X]AWV,(x) . (65) 

Here, there are L + 1 energy terms, Eo{x) and Vi{x) {i = 1,---,L), and A^^^ are the 
corresponding coupling constants for Ve{x). 

After integrating out the momentum degrees of freedom, the partition function of the 
system at fixed temperature T and parameters A is given by 

Z(T, X) = Jdx expi-(3E^ix)) = J dE^dV^ ■ ■ ■ dVL n{E^, V^i, ■ ■ ■ , Vl) exp {-^E^ , 

(66) 

where n{EQ, Vi, ■ ■ ■ , Vl) is the multidimensional density of states: 

n{Eo,Vi,---,VL) = J dx6{Eo{x)-Eo)6{Vi{x)-V)---6{VL{x)-VL) . (67) 

Here, the integration is replaced by a summation when x is discrete. 

The expression in Eq. fl65p is often used in simulations. For instance, in simulations of 
spin systems, Eq{x) and Vi{x) (here, L = 1 and x = {Si, 5*2, ■ ■ ■} stand for spins) can be 
respectively considered as the zero-field term and the magnetization term coupled with the 
external field A*^^^ (For Ising model, Eq = —JY.<i,j> SiSj, Vi = —Y.i Si, and A''^-' = h, i.e.. 
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external magnetic field.) In umbrella sampling [28] in molecular simulations, Eq{x) and 
Ve{x) can be taken as the original potential energy and the (biasing) umbrella potential 
energy, respectively, with the coupling parameter A*^^) (here, x = {q^, ■ ■ ■ , qjy} where q^. is 
the coordinate vector of the k-th particle and N is the total number of particles). For the 
molecular simulations in the isobaric-isothermal ensemble, Eo(x) and Vi{x) (here, L = 1) 
correspond respectively to the potential energy U and the volume V coupled with the 
pressure V. (Namely, we have x = {^i, ■ ■ ■ , ^at, V}, Eq = U, Vi = V, and A*^^-* = V, 
i.e., E^ is the enthalpy without the kinetic energy contributions.) For simulations in 
the grand canonical ensemble with particles, we have x = {qi, ■ ■ ■ , q^, TV}, and Eq{x) 
and Vi{x) (here, L = 1) correspond respectively to the potential energy U and the total 
number of particles N coupled with the chemical potential /i. (Namely, we have Eo = U, 
Vi = N, and A(i) = -^.) 

Moreover, going beyond the well-known ensembles discussed above, we can introduce 
any physical quantity of interest (or its function) as the additional potential energy term 
Ve. For instance, Ve can be an overlap with a reference configuration in spin glass systems, 
an end-to-end distance, a radius of gyration in molecular systems, etc. In such a case, we 
have to carefully choose the range of A^^'* values so that the new energy term X^^^Vg will 
have roughly the same order of magnitude as the original energy term Eq. We want to 
perform a simulation where a random walk not only in the Eq space but also in the Vg 
space is realized. As shown below, this can be done by performing a multidimensional 
REM, ST, or MUCA simulation. 

We first describe the multidimensional replica- exchange method (MREM) |101] . The 
crucial observation that led to this algorithm is: As long as we have M non-interacting 
replicas of the original system, the Hamiltonian H{q,p) of the system does not have to be 
identical among the replicas and it can depend on a parameter with different parameter 
values for different replicas. The system for the multidimensional REM consists of M non- 
interacting replicas of the original system in the "canonical ensemble" with M(= Mq x 
Ml X ■ ■ ■ X Ml) different parameter sets (m = 1, ■ ■ ■ , M), where A^ = (^mo' -^m) = 
(T^„, A« , ■ ■ ■ , A(,tJ) with mo = 1, ■ ■ ■ , Mo, = 1, • ■ ■ , (£ = 1, ■ ■ ■ , L). Because the 
replicas are non-interacting, the weight factor is given by the product of Boltzmann-like 
factors for each replica: 

Mo Ml Ml 

WWrem = n n ■ ■ ■ n exp (-^moEx^ ■ (68) 

mo=lmi=l mi^=l 

Without loss of generality we can order the parameters so that Ti < T2 < ■ ■ ■ < Tmo 
and Af^ < A2^'' < ■■■ < Aj^^, (for each I = 1,---,L). The multidimensional REM is 
realized by alternately performing the following two steps: 

1. For each replica, a "canonical" MC or MD simulation at the fixed parameter set is 
carried out simultaneously and independently for a certain steps. 

2. We exchange a pair of replicas i and j which are at the parameter sets A^ and 
Am+i, respectively. The transition probability for this replica exchange process is 
given by 

w{Am O A„+i) = min(l,exp(-A)) , (69) 

where we have 
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for T-exchange, and 



mo 



(71) 



for A^^-'-exchange (for one of i = 1, ■ ■ ■ , L). Here, g'*' and q^^^ stand for configuration 
variables for replicas i and j, respectively, before the replica exchange. 

We now consider the multidimensional simulated tempering (MST) which realizes a 
random walk both in temperature T and in parameters A [152^ 1153^ 1154] . The entire 
parameter set A = (T, A) = (T, A^^^, ■ ■ ■ , A^^)) become dynamical variables and both the 
configuration and the parameter set are updated during the simulation with a weight 
factor: 

W^mst(A) = exp (-(3Ex + /(A)) , (72) 

where the function /(A) = /(T, A) is chosen so that the probability distribution of A is 
fiat: 

Pmst(A) oc J dEodVi ■■■dVL n{Eo, ■ ■ ■ , Vl) exp {-PE^ + /(A)) = constant . (73) 

This means that /(A) is the dimensionless ( "Helmholtz" ) free energy: 

exp (-/(A)) = J dEodVi ■■■dVL n{Eo, 14, ■ ■ ■ , V^) expi-^E^) . (74) 

In the numerical work we discretize the parameter set A in M(= Mq x Mi x ■ ■ ■ x Ml) 
different values: A„, = (T^o,A^) = (T^^, A^ , . . . , A^^J), where mo = l,---,Mo,m£ = 
1, ■ ■ ■ , Ml {£ = 1, ■ ■ ■ , L). Without loss of generality we can order the parameters so that 
Ti < T2 < ■ ■ ■ < Tmo and A^^^ < < ■ ■ ■ < X^l]^ (for each i = 1, ■ ■ ■ , L). The free energy 

f{Am) is now written as fmo,rm,-,mL = /(^mo, A^^J, ■ ■ ■ , >^ml)- 

Once the initial configuration and the initial parameter set are chosen, the multidi- 
mensional ST is realized by alternately performing the following two steps: 



1. A "canonical" MC or MD simulation at the fixed parameter set A^ = (T^^, A^) = 
{Tmo ' ^ml ! ■ ■ ■ 5 ^ml ) is Carried out for a certain steps with the weight factor exp{—(3moEx 
(for fixed A^, /(A^) in Eq. (1721) is a constant and does not contribute). 

2. We update the parameter set A^ to a new parameter set A^ii in which one of 
the parameters in A^ is changed to a neighboring value with the configuration and 
the other parameters fixed. The transition probability of this parameter updating 
process is given by the following Metropolis criterion: 

w{Am Am±i) = min f 1, l^^^I^^^") = min (1, exp (-A)) . (75) 

\ l^MSTl-^mJ / 

Here, there are two possibilities for A^-i-i, namely, T-update and A*^^^-update. For 
T- update, we have A^ii = (^moiii ^m) with 

^ (/3mo±l l^mo) -^Am (/mo±l,mi,---,mi /mo,mi ,---,mi^ ) • C^^) 

For A^^^-update (for one of £ = 1, ■ ■ ■ , L), we have A^-i-i = (T^g, A^^^i) with 

where A ^^ =(■■■ A^^-^) X^'-K. A^^+i) ■ ■ and A =(■■■ A^^"^) A^ A(^+^) ■ ■ 
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We now describe the multidimensional multicanonical algorithm (MMUCA) which 
reahzes a random walk in the (L+l)-dimensional space of Eq[x) and Vi{x) {£ = 1, ■ ■ ■ , L). 

In the multidimensional MUCA ensemble, each state is weighted by the MUCA weight 
factor Wmmvca{Eo, Vi, ■ ■ ■ , Vl) so that a uniform energy distribution of Eq, Vi, ■ ■ ■, and 
Vl may be obtained: 



MMUCA 



{Eo, Vl,---, Vl) oc n{Eo, Vi,---, Vl)Wmmvca{Eo, Vi,---, Vl) = constant , (78) 



where u^Eq, Vi, - - - , Vl) is the multidimensional density of states. From this equation, we 
obtain 

1 

W^mmuca(^o, Vl,---, Vl) = exp (-/3a^MMUCA(^o, Vi,---, Vl)) oc 



n{Eo,Vi,---,VL) ' 

(79) 

where we have introduced an arbitrary reference temperature, = l/k-Q^a, and wrote 
the weight factor in the Boltzmann-like form. Here, the ^^multicanonical potential energij^ 
is defined by 

^mmuca(^o, Vl,---, Vl] T,) = ki,Ta lnn(Eo, Vi,---, Vl) . (80) 

The multidimensional MUCA MC simulation can be performed with the following Metropo- 
lis transition probability from state x with energy Ey^ = Eo + i;i=i A(^) V£ to state x' with 
energy E^ = E^' + Eti A^^V/ : 



w(x ^ X = mm 



' Wmmvca{Eo',Vi',---,Vl') \ ^ . ( n{Eo,Vi,---,VL) 
, ' H^MMUCA(i^o, Vl,---, Vl) ) ^ ' niEo', V, Vl') 



A MD algorithm in the multidimensional MUCA ensemble also naturally follows from 
Eq. (179|) . in which a regular constant temperature MD simulation (with T = T^) is 
performed by replacing the total potential energy E^ by the multicanonical potential 
energy -Emmuca in Eq. (Il2]): 



dEMMVCA{Eo,Vi, - - - ,VL;Ta) s 

= d^, sp^- ^''^ 

We remark that the random walk in Eq and in Ve for the MUCA simulation corresponds 
to that in (3 and in /JA'-^^ for the ST simulation: 

/ Eq i — > /3 , , , 

j V ^/3AW, {i=l,---,L). 

They are in conjugate relation. 



3.2 Weight Factor Determinations for Multidimensional ST and 
MUCA 

Among the three multidimensional generalized-ensemble algorithms described above, only 
MREM can be performed without much preparation because the weight factor for MREM 
is just a product of regular Boltzmann-like factors. On the other hand, we do not know 
the MST and MMUCA weight factors a priori and need to estimate them. As a simple 
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method for these weight factor determinations, we can generahze the REST and REMUCA 
presented in the previous subsections to muhidimensions. 

Suppose we have made a single run of a short MREM simulation with M(= Mq x Mi x 
■ ■ ■ X Ml) replicas that correspond to M different parameter sets {fri = 1, ■ ■ ■ , M). 
Let Amo,mi,-,mi(-£^o, Vl,- - ■ , Vl) and nmo,mu-,mL be respectively the (L + l)-dimensional 
potential-energy histogram and the total number of samples obtained for the m-th pa- 
rameter set Am = {Tmo, ■ ■ ■ , A^]). The generalized WHAM equations are then given 
by 

E Nmo,mi,-,mL{Eo',Vi, - - - ,Vl) 

n{Eo, V,---, Vl) = ^ _ , (84) 

mo,mi,---,mL 

and 

exp{-fmo,mu-,mJ = E n{Eo, Vl, - - - , Vl) exp ^rnoE x^) ■ (85) 

Eo,Vi,-,Vl 

The density of states n{Eo, Vi, - - - , Vl) (which is inversely proportional to the MMUCA 
weight factor) and the dimensionless free energy fmo,mi,- -,mL (which is the MST parameter) 
are obtained by solving Eqs. fl84p and (153]) self-consistently by iteration. 



3.3 Expectation Values of Physical Quantities 

We now present the equations to calculate ensemble averages of physical quantities with 
any temperature T and any parameter A values. 

After a long production run of MREM and MST simulations, the canonical expectation 
value of a physical quantity A with the parameter values A^ {m = 1, - - - , M), where A^ = 
(Tmo, Xm) = (Tmo, , ■ ■ ■ , A^]) with mo = 1, ■ ■ ■ , Mq, me = 1, - - - , Mg {£ = 1, - - - , L), and 
M(= Mo X Ml X ■ ■ ■ X Ml), can be calculated by the usual arithmetic mean: 

^ rim 

<A>^ X = —Y.A{xm{k)) , (86) 

where Xm{k) {k = 1, - - - the configurations obtained with the parameter values 

Am {fn = 1,---,M), and rim is the total number of measurements made with these 
parameter values. The expectation values of A at any intermediate T (= I/Zcb/?) and any 
A can also be obtained from 

^ A{Eo, Vl,---, VL)n{Eo, Vi,---, Vl) exp [-/3Ex 

E n{Eo,Vi,---,VL)exp{-^Ex) 

Eoyi,-yL 

where the density of states n{EQ,Vi, - - - ,Vl) is obtained from the multiple-histogram 
reweighting techniques. Namely, from the MREM or MST simulation, we first obtain 
the histogram A^mo,mi,- -,mi(-E'o! Vi, - - - , Vl) and the total number of samples ?^mo,mi,---,mi 
in Eq. flM]) . The density of states n{EQ,Vi, - - - ,Vl) and the dimensionless free energy 
fmo,mi,---,mL then obtained by solving Eqs. flM|) and f lSS]) self-consistently by itera- 
tion. Substituting the obtained density of states n{Eo, Vi, - - - , Vl) into Eq. ( 187|) . one can 
calculate the ensemble average of the physical quantity A at any T and any A. 
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Moreover, the ensemble average of the physical quantity A (including those that cannot 
be expressed as functions of Eq and Vj (£ = 1, ■ ■ ■ , L) ) can be obtained from the "trajec- 
tory" of configurations of the production run |153j . Namely, we first obtain fmo,mi,---,mL 
for each (mo = 1, ■ ■ ■ , Mq, mi = 1, ■ ■ • , Mi, ■ ■ ■ , m^ = 1, ■ ■ ■ , Ml) by solving Eqs. (184|) and 
(185|) self-consistently, and then we have 



Mo Ml 

E ■■■ E EM^r 



exp (-/3Exi 



Mo Ml 

mo=l m_t = l x'„ ^ ^ _ „„„ , ^- w ^ ,™ 



E ' ' ' E ^no,---,nL 6Xp (/no, ^"0-^A„ 



^ _ _ exp {-PE^jxr, 

2^ ' ' ' ^ ^ Afo AfL 

E E ^no,---,nL exp (/ng^...^„^ — /^np-E^^ (Xn 

"0=1 »T.L=1 

(88) 

where the configurations obtained at Am = (Tm(j,Am) = (T^q, A^|, ■ ■ ■ , A^]). 

Here, the trajectories are stored for each separately. 

For the MMUCA simulation with the weight factor Wmmuca(-E'o, ■ ■ ■ , Vl), the expec- 
tation values of A at any T (= 1//cb/3) and any A can also be obtained from Eq. ( 187|) by 
the single-histogram reweighting techniques as follows. Let A'mmuca(-E'o, ^i, ■ ■ ■ , ^l) be 
the histogram of the distribution of Eq,Vi, ■ ■ ■ ,Vl, Pmmuca(-E'o) Vi, ■ ■ ■ ,Vl), obtained by 
the production run. The best estimate of the density of states u^Eq, Vi, ■ ■ ■ , Vl) is then 
given by 

n{EQ Vl ■ ■ ■ Vl) = Vi, - ■ ■ ,Vl) ^gg^ 

M^mmuca(-E'o5 ■ ■ ■ 5 Vl) 

Moreover, the ensemble average of the physical quantity A (including those that cannot 
be expressed as a function of Eq and Vg {i = 1, ■ ■ ■ , L)) can be obtained as long as one 
stores the "trajectory" of configurations Xk from the production run. We have 



J2 Mxk)W^li^cAiEo{xk), ■ ■ ■ , VL{xk)) exp (^-l3E)^{xk) 

< A >^ x= ^ • (90) 

E ^MMUcA(^o(a;fc), ■ ■ ■ , VL{xk)) exp (^-^E^ixk) 



k=l 



Here, Xk is the configuration at the k-th MC (or MD) step and Ug is the total number of 
configurations stored. 



3.4 Multidimensional Generalized-Ensemble Algorithms for the 
Isobaric-Isothermal Ensemble 

As examples of the multidimensional formulations in the previous subsections, we present 
the generalized-ensemble algorithms for the isobaric-isothermal ensemble (or, the NPT 
ensemble) |155j . Let us consider a physical system that consists of N atoms and that 
is in a box of a finite volume V. The states of the system are specified by coordinates 
q = {q^, q2, - -, q^} momenta p = {Pi,P2, ■ ■ ■ ,Pn} of the atoms and volume V of 
the box. The potential energy E(q, V) for the system is a function of q and V. 
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In the isobaric-isothermal ensemble |157[ I158[ I163[ 1170] the probability distribution 

Pj<{pt{E, V; T, V) for potential energy E and volume V at temperature T and pressure V 
is given by 



Pnpt(^, V; T, V) oc niE, V)W^pt{E, V; T, V) = n{E, V)e-f^ 



(91) 



Here, the density of states n{E, V) is given as a function of both E and V, and "H is the 
"enthalpy" (without the kinetic energy contributions): 



n = E + vv 



(92) 



This weight factor produces an isobaric-isothermal ensemble at constant temperature (T) 
and constant pressure (V). Note that this is a special case of the general formulations in 
Eq. ([65]) with L = l, Eo = E,Vi = V, and A^^) = V. 

In order to perform the isobaric-isothermal MC simulation |170] . we perform Metropo- 
lis samphng on the scaled coordinates a = {cti, a2, ■ ■ ■ ,(Tn} where cr^ = V'^^^q^ (k = 
1, 2, ■ ■ ■ , A^) (q^j are the real coordinates) and the volume V (here, the particles are placed 
in a cubic box of a side of size V~^^^). The trial moves from state x with the scaled 
coordinates a with volume V to state x' with the scaled coordinate a' and volume V 
are generated by uniform random numbers. The enthalpy is accordingly changed from 
T-L^E^a, V), V) to H'^E^a', V), V) by these trial moves. The trial moves will be accepted 
with the following Metropolis criterion: 



w{x x') = min (l,exp[-/3{H' - H - iVA;BTln(V7V)}]) 



(93) 



where N is the total number of atoms in the system. 

As for the MD method in this ensemble, we just present the Nose-Andersen algo- 
rithm |157[ I158[ 1163] . The equations of motion in Eqs. ( lTT]) - (fT4|) are now generalized as 
follows: 
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where Ai is the artificial mass associated with the volume, Py is the conjugate momentum 
for the volume, and the "instantaneous pressure" V{t) is defined by 



(102) 



In REM simulations for the NPT ensemble, we prepare a system that consists of 
Mt X M-p non-interacting replicas of the original system, where My and Mp are the 
number of temperature and pressure values used in the simulation, respectively. The 
replicas are specified by labels i [i = 1,2, Mt x M-p), temperature by mo (mo = 
1, 2, ■ ■ ■ , Mt) and pressure by mi (mi = 1, 2, ■ ■ ■ , Mp). 

To perform REM simulations, we carry out the following two steps alternately: (1) per- 
form a usual constant NPT MC or MD simulation in each replica at assigned temperature 
and pressure and (2) try to exchange the replicas. If the temperature (specified by mo and 
no) or pressure (specified by mi and ni) between the replicas is exchanged in Step 2, the 
transition probability from X = {■ ■ ■ , (ciW, V^; T^^, P^J, ■ ■ ■ , (a^, Vl^l; T„„, ■ ■ ■} to 
X' = {■■■, (ctH, VW; Tn„ PnJ, ■ ■ ■ , (atJl, V^^\ T^,,V^,), ■ ■ ■} at the trial is given by [mi [II] 

wrem(^ ^0 = min [1, exp(-AREM)] , (103) 

where 



AREM = (/3mo-/3no) ^(a[^lV[^l)-E(0-W,V'*') +(/3™o^^;ni-/3noPnJ(V^^'-V'*'). (104) 



In ST simulations for the NPT ensemble, we introduce a function f{T,V) and use 
a weight factor Wst{E,V;T,V) = exp[-/3(E + VV) + f{T,V)] so that the distribution 
function P^t^{T, V) of T and V may be uniform: 

POO p 

Pst(T, V)(x dV dq WsT[E{q, V), V; T, V] = constant . (105) 
Jo Jv 

From Eq. fllOSp . it is found that f{T,V) is formally given by 

/(T, P) = - In 1 dV 1^ dq exp [-(3 {E{q, V) + PV)] } , (106) 

and this function is the dimensionless Gibbs free energy except for a constant. 

To perform ST simulations, we again discretize temperature and pressure into Mo x Mi 
set of values (T^p, Pmi) ("^o = 1? " " " ? Mq, mi = 1, ■ ■ ■ , Mi). We carry out the following 
two steps alternately: (1) perform a usual constant NPT MC or MD simulation and (2) 
try to update the temperature or pressure. In Step 2 the transition probability from the 
state X = {a, V; Tmo,Vmi} to the state X' = {a, V; Tno,Vni} is given by 

wsTiX ^ X') = min [1, exp(-AsT)] , (107) 

where 

AsT = (/3no - l3mo)E{a, V) + ilSnoVn, " /^mo^mJV - (/„„,„, - fm,,,mj • (108) 

We remark that when we perform MD simulations with REM and ST, the momenta 
should be rescaled if the replicas are exchanged for the temperature in REM and the 
temperature is updated in ST as shown above in the previous subsections. 
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From the production run of REM or ST simulations in the NPT ensemble, we can 
calculate isobaric-isothermal averages of a physical quantity A at iTmo,'PmJ {mQ = 
1, ■ ■ ■ , Mq, mi = 1, ■ ■ ■ , Ml) by the usual arithmetic mean: 

rim 

<A>Tm,,Vm=—J2^Mk)) , (109) 

where Xm{k) {k = 1, - ■ ■ ,nm) are the configurations obtained with the parameter values 
[Tmo , Vmi ) and rim is the total number of measurements made with these parameter 
values. The expectation values of A at any intermediate temperature T (= l/k^P) and 
any intermediate pressure V can also be obtained from 

J2 ME, V)n{E, V) exp {-I3{E + VV)) 
<A>TP = ^'^^ , (110) 

E,V 

where the density of states n{E, V) is obtained from the multiple-histogram reweighting 
techniques. Namely, from the REM or ST simulation, we first obtain the histogram 
Nmo,mi{E, V) and the total number of samples rimo^mi- The density of states n{E, V) and 
the dimensionless free energy fmo,mi are then obtained by solving the following equations 
self-consistently by iteration (see Eqs. ( 184|) and ( 185|) above): 

niE, V) = "-^ , (111) 

/ , ^mo,mi GXp (/"mo, mi f^moi^E -\- T^mi^)) 
mo, mi 

and 

exp(-/^o,™J = ^n(E,V)exp(-/3^„(E + P^,V)) . (112) 

E,V 

Substituting the obtained density of states n{E,V) into Eq. flllOp . one can calculate the 
ensemble average of the physical quantity A at any T and any V. 

We now introduce the multicanonical algorithm into the isobaric-isothermal ensemble 
and refer to this generalized-ensemble algorithm as the multiharic-multithermal algorithm 
(MUBATH) [70]-[73]. The molecular simulations in this generalized ensemble perform 
random walks both in the potential energy space and in the volume space. 

In the MUBATH ensemble, each state is sampled by the MUBATH weight factor 
Wmbt(-E', V) = exp{— /3a'Hmbt(-E, V)} (Hmbt IS referred to as the multibaric-multithermal 
enthalpy) so that a uniform distribution in both potential energy E and volume V is 
obtained 1701: 



P„ibt(^, V) cx n{E, V)iyinbt(^, V) = n{E, V) exp{-/3„-H^bt(^, V)} = constant , (113) 

where we have chosen an arbitrary reference temperature, Tq = l/k-B^a- 

The MUBATH MC simulation can be performed by replacing "H by "Hmbt in Eq- (|93l) : 

w{x ^ x') = min (1, exp[-/3,{-HUt - ^mbt - Nk,,Ta ln(V7V)}]) , (114) 

In order to perform the MUBATH MD simulation, we just solve the above equations 
of motion (Eqs. (I94l) - fll01l) ) for the regular isobaric-isothermal ensemble (with arbitrary 
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reference temperature T = To), where the enthalpy H is replaced by the multibaric- 
multithermal enthalpy "Hmbt in Eqs. f l^ and f llOOp [72] . 

In order to calculate the isobaric-isothermal-ensemble averages, we employ the single- 
histogram reweighting techniques [17] . The expectation value of a physical quantity A at 
any T and any V is obtained by substituting the following density of states into Eq. f lllOp : 

where iVmbt(-E, V) is the histogram of the probability distribution P^yyt{E, V) of potential 
energy and volume that was obtained by the MUBATH production run. 



4 EXAMPLES OF SIMULATION RESULTS 

We tested the effectiveness of the generalized-ensemble algorithms by using a system of 
a 17- residue fragment of ribonuclease Ti |17H I1U8] . It is known by experiments that 
this peptide fragment forms a- helical conformations |171] . We have performed a two- 
dimensional REM simulation and a two-dimensional ST simulation. In these simulations, 
we used the following energy function: 

Ex = Eo + XEsoL, (116) 

where we set L = l,Vi = Esql, and A^^^ = A in Eq. fl65|) . Here, Eq is the potential energy 
of the solute and -Esol is the solvation free energy. The parameters in the conformational 
energy as well as the molecular geometry were taken from ECEPP/2 [172^ 1173^ 1174] . 

The solvation term Esol is given by the sum of terms that are proportional to the 
solvent-accessible surface area of heavy atoms of the solute [175] . For the calculations of 
solvent-accessible surface area, we used the computer code NSOL |176] . 

The computer code KONF90 [71|S] was modified in order to accommodate the generalized- 
ensemble algorithms. The simulations were started from randomly generated conforma- 
tions. We prepared eight temperatures (Mq =8) which are distributed exponentially 
between Ti = 300 K and Tmq = 700 K (i.e., 300.00, 338.60, 382.17, 431.36, 486.85, 549.49, 
620.20, and 700.00 K) and four equally-spaced A values (Mi = 4) ranging from to 1 
(i.e., Ai = 0, A2 = 1/3, A3 = 2/3, and A4 = 1) in the two-dimensional REM simulation 
and the two-dimensional ST simulation. Simulations with A = (i.e.. Ex = Eq) and with 
A = 1 (i.e.. Ex = Eq + -Esol) correspond to those in gas phase and in aqueous solution, 
respectively. 

We first present the results of the two-dimensional REM simulation. We used 32 
replicas with the eight temperature values and the four A values given above. Before 
taking the data, we made the two-dimensional REM simulation of 100000 MC sweeps 
with each replica for thermalization. We then performed the two-dimensional REM sim- 
ulation of 1000000 MC sweeps for each replica to determine the weight factor for the 
two-dimensional ST simulation. At every 20 MC sweeps, either T-exchange or A-exchange 
was tried (the choice of T or A was made randomly). In each case, either set of pairs of 
replicas ((1,2),...,(M— 1,M)) or ((2,3),...,(M,1)) was also chosen randomly, where M is 
Mq and Mi for T-exchange and A-exchange, respectively. 

In Fig. 1 we show the time series of labels of T^o (i.e., mo) and Ami (i-^-, ^i) for one of 
the replicas. The replica realized a random walk not only in temperature space but also 
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(a) 




Figure 1: Time series of the labels of T^o, mo, (a) and Xm.-^, mi, (b) as functions of MC 
sweeps, and that of both mo and mi for the region from 400000 MC sweeps to 700000 
MC sweeps (c). The results were from one of the replicas (Replica 1). In (a) and (b), 
MC sweeps start at 100000 and end at 1100000 because the first 100000 sweeps have been 
removed from the consideration for thermalization purpose. 



in A space. The behavior of T and A for other replicas was also similar (see Ref. [154] ). 
From Fig. 1, one finds that the A-random walk is more frequent than the T-random walk. 




200000 600000 1e+06 200000 600000 1e+06 

MC sweeps MC sweeps 



Figure 2: Time series of the temperature T (a), total energy -Etot (b), conformational 
energy Eq (c), solvation free energy -Esol (d), and end-to-end distance D (e) for the 
same replica as in Fig. 1. The temperature is in K, the energy is in kcal/mol, and the 
end-to-end distance is in A. 

We also show the time series of temperature T, total energy -Etot, conformational 
energy Eq, solvation free energy Esol, and end-to-end distance D for the same replica 
in Fig. 2. From Figs. 2(a) and 2(e), we find that at lower temperatures the end-to-end 
distance is about 8 A, which is the length of a fully a-helical conformation and that at 
higher temperatures it fluctuates much for a range from 7 A to 14 A. It suggests that 
a-helix structures exist at low temperatures and random-coil structures occur at high 
temperatures. There are transitions from/to a-helix structures to/from random coils 
during the simulation. It indicates that the REM simulation avoided getting trapped in 
local-minimum-energy states and sampled a wide conformational space. 

The canonical probability distributions of -Etot and -Esol at the 32 conditions ob- 
tained from the two-dimensional REM simulation are shown in Fig. 3. For an optimal 
performance of the REM simulation, there should be enough overlaps between all pairs 
of neighboring distributions, which will lead to sufficiently uniform and large acceptance 
ratios of replica exchanges. There are indeed ample overlaps between the neighboring 
distributions in Fig. 3. 
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Figure 3: Contour curves and histograms of distributions of the total energy Etqt and 
the solvation free energy Esoh ((a) and (b)) from the two-dimensional REM simulation. 



We now use the results of the two-dimensional REM simulation to determine the 
weight factors for the two-dimensional ST simulation by the multiple-histogram reweight- 
ing techniques. Namely, by solving the generalized WHAM equations in Eqs. f lH^ and 
f l85|) with the obtained histograms at the 32 conditions (see Fig. 3), we obtained 32 values 
of the ST parameters fmo,m.i {rno = 1, ■ ■ ■ , 8; mi = 1, ■ ■ ■ , 4). 

After obtaining the ST weight factor, Wst 
carried out the two-dimensional ST simulation of 1000000 MC sweeps for data collection 
after 100000 MC sweeps for thermalization. At every 20 MC sweeps, either or 
was respectively updated to T^oii or Ami±i (the choice of T or A update and the choice 
of ±1 were made randomly). 

We show the average total energy, average conformational energy, average A x Esol, 
and average end-to-end distance in Fig. 4. The results are in good agreement with those 
of the REM simulation (data not shown). 




"300 400 500 600 700 300 400 500 600 700 



(a) T 




300 400 500 600 700 300 400 500 600 700 



(c) fd) T 

Figure 4: The average total energy (a), average conformational energy (b), average of 
A X Esoh (c), and average end-to-end distance (d) with all the A values as functions of 
temperature. The lines colored in red, green, blue, and purple are for Ai, A2, A3, and A4, 
respectively. They are in order from above to below in (a) and (c) and from below to 
above in (b) and (d). 



We found that the results of the two-dimensional ST simulation are in complete agree- 
ment with those of the two-dimensional REM simulation for the average quantities. The 
only difference between the two simulations is the number of replicas. In the present 
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simulation, while the REM simulation used 32 replicas, the ST simulation used only one 
replica. Hence, we can save much computer power with ST. 

A second example of our multidimensional generalized-ensemble simulations is a pres- 
sure ST (PST) simulation in the isobaric-isothermal ensemble |155] . This simulation 
performs a random walk in one-dimensional pressure space. The system that we simu- 
lated is ubiquitin in explicit water. This system has been studied by high pressure NMR 
experiments and known to undergo high-pressure denaturations \177\ 1178] . Ubiquitin has 
76 amino acids and it was placed in a cubic box of 6232 water molecules. Temperature 
was fixed to be 300 K throughout the simulations, and we prepared 100 values of pressure 
ranging from 1 bar to 10000 bar. Temperature and pressure were controlled by Hoover- 
Langevin method |179] and particle mesh Ewald method |180[ 1181] were employed for 
electrostatic interactions. The time step was 2.0 fsec. The force field CHARMM22 |182] 
with CMAP MM and TIP3P water model [IHSl [IH2] were used, and the program 
package NAMD version 2.7b3 |186j was modified to incorporate the PST algorithm. 

We first performed 100 independent conventional isobaric-isothermal simulations of 4 
nsec with T = 300 K (i.e., Mq = 1) and 100 values of pressure (i.e.. Mi = 100). Using 
the obtained histogram Nmo,mi{E,V) of potential energy and volume distribution, we 
obtained the ST parameters /mo, mi by solving the WHAM equations in Eqs. f illip and 
f lll2p . We then performed the PST production of 500 nsec and repeated it 10 times with 
different seeds for random numbers (so, the total simulation time for the production run 
is 5.0 /isec). 

In Fig. 5 we show the time series of pressure and potential energy during the PST 
production run. 




Time (ns) Time (ns) 

Figure 5: Time series of pressure (left) and potential energy (right) during the PST 
production run. 

In the Figure we see a random walk in pressure between 1 bar and 10000 bar. A 
random walk in potential energy is also observed and it is anti-correlated with that of 
pressure, as it should be. 

We calculated the fluctuations < d"^ > — < d of the distance d between pairs of 
C° atoms. The results are shown in Fig. 6. 

We see that large fluctuations are observed between residues around 7-10 and around 
20-40, which are in accord with the experimental results |177[ 1178] . 

The fiucutating distance corresponds to that between the turn region of the /3-hairpin 
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Figure 6: Fluctuations of distance between pairs of atoms that was calculated from 
the PST production run. 



and the end of the a-helix as depicted in Fig. 7. While at low pressure this distance is 
small, at high pressure it is larger and water comes into the created open region. 



5 CONCLUSIONS 

In this article we first introduced three well-known generalized-ensemble algorithms, 
namely, REM, ST, and MUCA, which can greatly enhance conformational sampling of 
biomolecular systems. We then presented various extensions of these algorithms. Ex- 
amples are the general formulations of the multidimensional REM, ST, and MUCA. We 
generalized the original potential energy function Eq by adding any physical quantities Vi 
of interest as a new energy term with a coupling constant A*^^^(^ = 1, ■ ■ ■ , L). The simu- 
lations in multidimensional REM and multidimensional ST algorithms realize a random 
walk in temperature and = 1, ■ ■ ■ , L) spaces. On the other hand, the simulation in 

multidimensional MUCA algorithms realizes a random walk in Eq,Vi, ■ ■ ■ ,Vl spaces. 

While the multidimensional REM simulation can be easily performed because no 
weight factor determination is necessary, the required number of replicas can be quite 
large and computationally demanding. We thus prefer to use the multidimensional ST or 
MUCA, where only a single replica is simulated, instead of REM. However, it is very dif- 
ficult to obtain optimal weight factors for the multidimensional ST and MUCA. Here, we 
have proposed a powerful method to determine these weight factors. Namely, we first per- 
form a short multidimensional REM simulation and use the multiple-histogram reweight- 
ing techniques to determine the weight factors for multidimensional ST and MUCA sim- 
ulations. 

The multidimensional generalized-ensemble algorithms that were presented in the 
present article will be very useful for Monte Carlo and molecular dynamics simulations of 
complex systems such as spin glass, polymer, and biomolecular systems. 
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Figure 7: Snapshots of ubiquitin during the PST production run at low pressure (left) 
and at high pressure (right). 
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